MBD Estadística - Practica 1

A continuación se expone el código y comentarios resultados de la realización de la práctica 1.

0. Imports

Importar librerias

library(tidyverse)
library(tidyr)
library(dplyr)
library(knitr)
library(ggplot2)
library(stringr)
library(reshape)
library(minerva)
library(heatmaply)
library(kableExtra)
library(factoextra)
library(ggbiplot)
library(randomForest)
library(ROCR)
library(ggcorrplot)
library(ModelMetrics)

Importar datos

Aunque se nos entrega el fichero de train entero y podriamos obtener la variable respuesta, se trataran cómo si de un caso real se tratara, y todo el EDA y el training se usara solo el p1_train. El p1_test se usara únicamente para validar los resultados al final de todo de la practica.

train <- read.csv("./data/train.csv", head = TRUE, stringsAsFactors=TRUE)

p1_train <- read.csv("./data/p1_train.csv", head = TRUE, sep=";", stringsAsFactors=TRUE)

p1_test <- read.csv("./data/p1_test.csv", sep=";", head = TRUE, stringsAsFactors=TRUE)
p1_test_fin <- read.csv("./data/p1_test.csv", sep=";", head = TRUE, stringsAsFactors=TRUE)

El df p1_test_fin se crea para mantener los datos base sin transformaciones para la generación final del resultado.

1. Preprocesado de los datos

En este apartado se va hacer el análisi exploratorio de los datos.

head(p1_train)

Eliminamos la variable id del dataset ya que es el indice, y no aporta información.

p1_train <- p1_train %>% select(-id)

Summary de los datos

summary(p1_train)

      year           hour           season         holiday        workingday        weather          temp           atemp      
 Min.   :2011   Min.   : 0.00   Min.   :1.000   Min.   :0.000   Min.   :0.0000   Min.   :1.00   Min.   : 0.82   Min.   : 0.76  
 1st Qu.:2011   1st Qu.: 6.00   1st Qu.:2.000   1st Qu.:0.000   1st Qu.:0.0000   1st Qu.:1.00   1st Qu.:13.94   1st Qu.:16.66  
 Median :2011   Median :12.00   Median :3.000   Median :0.000   Median :1.0000   Median :1.00   Median :20.50   Median :24.24  
 Mean   :2011   Mean   :11.57   Mean   :2.506   Mean   :0.029   Mean   :0.6773   Mean   :1.41   Mean   :20.27   Mean   :23.70  
 3rd Qu.:2012   3rd Qu.:18.00   3rd Qu.:4.000   3rd Qu.:0.000   3rd Qu.:1.0000   3rd Qu.:2.00   3rd Qu.:26.24   3rd Qu.:31.06  
 Max.   :2012   Max.   :23.00   Max.   :4.000   Max.   :1.000   Max.   :1.0000   Max.   :3.00   Max.   :41.00   Max.   :45.45  
    humidity        windspeed          count      
 Min.   :  0.00   Min.   : 0.000   Min.   :  1.0  
 1st Qu.: 46.00   1st Qu.: 7.002   1st Qu.: 41.0  
 Median : 62.00   Median :12.998   Median :145.0  
 Mean   : 61.77   Mean   :12.802   Mean   :191.4  
 3rd Qu.: 77.00   3rd Qu.:16.998   3rd Qu.:283.0  
 Max.   :100.00   Max.   :56.997   Max.   :977.0  

1.0 Nan handling

Se puede apreciar que no encontramos ningún nan en las columnas del dataset, haciendo que no tengamos que tratar con ellos.

p1_train %>% group_by() %>% summarise_all(funs(sum(is.na(.))))

1.1 Duplicates handling

Se eliminan las filas duplicadas del dataset.

p1_train <- p1_train %>%
                distinct()
p1_train

1.2. Separar datos entre numericos y categoricos

Para las variables que son categoricas, aunque su representación sea numerica debemos tratarlas cómo categoricas a nivel dato para obtener el resultado deseado.

Las variables categoricas del dataset segun el enunciado son: workingday, holiday, weather y season.

p1_train <- p1_train %>% mutate(workingday = as.factor(workingday),
                                   holiday = as.factor(holiday),
                                   weather = as.factor(weather),
                                   season = as.factor(season))

2. Exploratory Data Analysis

Eliminamos la variable target del dataset para hacer un estudio de las columnas predictivas. Posteriormente se hará el estudio de la variable target independientemente.

predictors <- p1_train %>% select(-count)

2.0. Hacer boxplots de los campos y mirar outliers en los datos

Estudio de las variables predictivas numericas.

meltPred <- predictors %>% select(where(is.numeric)) %>% melt()

Using  as id variables

meltPred %>%
 ggplot(aes(factor(variable), value)) +
   geom_violin(width=1, color = "gray", alpha = 0.2, fill = 'green' ) +
    geom_boxplot(width=0.3, color="black", fill='blue', alpha=0.3, outlier.colour="red", outlier.shape=8,
             outlier.size=1, notch=TRUE) + facet_wrap(~variable, scale="free")

Notch went outside hinges
ℹ Do you want `notch = FALSE`?

Estos gráficos representan un boxplot y un violinplot sobre todas las variables numericas del dataset. Podemos apreciar:

• La columna year contiene únicamente dos valores 2011 y 2012.

• La columna hour se ditribuye entre el 0 y el 24, con una media situada en los 12. No encontramos outliers.

• La columna temp y la columna atemp siguen distribuciones muy similares, la media esta alrededor del 25. No encontramos outliers.

• La columna humidity se distribuye entre el 0 y el 100 con una media alrededor de los 60. No encontramos otuliers.

• La columna windspeed se distribuye entre el 0 y el 60 con una media alrededor de los 15. Encontramos numersos outliers cuando el valores es mayor de 30.

2.1 Matriz de correlaciones de pearson de las variables numericas

corrp <- p1_train %>% select(where(is.numeric)) %>%  cor(method = 'pearson')
corrp %>% ggcorrplot(hc.order = TRUE, type = 'lower', lab=TRUE,
                      outline.col = "white",
                       ggtheme = ggplot2::theme_gray,
                       colors = c("#6D9EC1", "white", "#E46726")
                      )

Se puede apreciar cómo la variable atemp y temp estan muy correlacionadas, de ambas, seleccionaremos solamente la columna temp y eliminaremos la columna atemp para una mejora en la prediccion. Con la variable target count, estan directamente relacionadas las columnas temp, year, hour y windspeed y la columna humidity esta indirectamente realcionada.

Eliminamos la columna atemp del dataset.

p1_train <- p1_train %>% select(-atemp)
p1_test <- p1_test %>% select(-atemp)
predictors <- predictors %>% select(-atemp)

2.2 matriz de correlaciones lineales y no lineales con la variable target ‘count’

Calculamos la correlación no lineal MIC de las columnas con la variable target.

numpred <- predictors %>% select(where(is.numeric)) 
corrs <- data.frame(MIC = mine(numpred, y = p1_train$count, alpha = 0.7)) %>% select(MIC.Y)
corrs$Pearson <- cor(numpred, y = p1_train$count, method = 'pearson')
rownames(corrs) <- rownames(corrs$Pearson)
corrs

Podremos apreciar que las correlaciones no linealas calculadas son (en proporcion) similares a las correlaciones de Pearson, entendiendo entonces que no hay correlaciones obviadas en el grafico de correlaciones.

cbind(corrs$MIC.Y, corrs$MIC.Y ) %>% heatmaply_cor(xlab = "",
              ylab = "Columnas", main = "Correlacion MIC de las columnas numericas con la variable target 'count'", cellnote = cbind(corrs$MIC.Y, corrs$MIC.Y ), k_col = 1, k_row = 1)

NA

El gràfico superior tiene 2 columnas pero se debe considerar una sola (estan repetidas), ha sido la unica manera que hemos encontrado de hacer el gráfico del MIC.

2.3 PCA

Realizamos un PCA sobre los datos para extraer los componentes principales y estudiarlos.

pca <- predictors %>% select(where(is.numeric)) %>% sample_n(200) %>%
  prcomp(scale. = TRUE, center = TRUE) #Se debe hacer un sample para posteriormente poder apreciar graficamente la distribución PCA.
summary(pca)

Importance of components:
                         PC1    PC2    PC3    PC4    PC5
Standard deviation     1.257 1.1315 0.9295 0.8425 0.7523
Proportion of Variance 0.316 0.2560 0.1728 0.1419 0.1132
Cumulative Proportion  0.316 0.5721 0.7449 0.8868 1.0000

pca %>% 
  ggbiplot::ggbiplot(scale = 1)

Gràfico que representa la ponderación de las variables utilizadas para el PCA.

3. Analisis exploratorio de la variable target ‘count’

La variable target ‘count’ es una variable numerica que representa el numero de bicicletas alquiladas en las grandes ciudades. La variable tiene una media de 190 y un std de 182, un valor bastante elevado que indica mucha variación en la variable.

corrs$Pearson

                [,1]
year       0.2623084
hour       0.4037404
temp       0.3971578
humidity  -0.3245222
windspeed  0.1078587

Correlaciones de pearson de las columnas predictivas con la variable ‘count’.

mean(p1_train$count)

[1] 191.5154

sd(p1_train$count)

[1] 182.139

3.1 Distribución

El gráfico siguiente es un boxplot y un violin plot sobre la variable target ‘count’. Por el violin plot se aprecia que la mayoria de valores se encuentran cercanos al 0 pero por el boxplot se visualiza que los valores altos del dataset hacen que la media suba considerablemente. En el boxplot podemos apreciar también diversos outliers que sobresalen a partir del valor 550.

p1_train %>%
  ggplot( aes(x=count, y=count)) +
    geom_violin(width=1, color = "gray", alpha = 0.2, fill = 'green' ) +
    geom_boxplot(width=0.3, color="black", fill='blue', alpha=0.3, outlier.colour="red", outlier.shape=8,
             outlier.size=1, notch=TRUE) +
    theme(
      legend.position="none",
      plot.title = element_text(size=11)
    ) +
    ggtitle("Violin y box plot de los valores de la variable target ") +
    xlab("") + theme_classic() 

3.2 Histograma y densidad

Partiendo de la base que la variable target es continua, podemos apreciar que la frequencia de la misma sigue una distribución parecida a la exponencial, ya que cómo más se aleja del número 0 menos frequencia encontramos. Habrá que tener en cuenta entonces, que el modelo va a tender a predecir valores bajos en la mayoria de los casos.

p1_train %>%
ggplot( aes(x=count)) + 
 geom_histogram(aes(y=..density..), colour="black", fill="white", binwidth = 20)+
 geom_density(alpha=.2, fill="#FF6666") + theme_classic()

4. Prediccion

4.1 Normalización y preparación de los datos

Hacemos un escalado sobre los datos para normalizar.

p1_train <- p1_train %>% mutate(across(where(is.numeric), scale))

Separamos el dataset para train y test, para poder hacer validaciones sobre el modelo.

p <- 0.8
n <- nrow(p1_train)
set.seed(12345)
train.sel <- sample(c(FALSE, TRUE), n, rep = TRUE, prob=(c(1-p,p)))
train <- p1_train[train.sel,]
test <- p1_train[!train.sel,]

train

4.2 Linear model

Entrenamos un modelo lineal con los datos del train

mod_lm <- lm(count~., train) 

Calculamos el mse sobre el train

mean(mod_lm$residuals^2) #mse

[1] 20124.92

Calculamos el msle sobre el train

mean(log(mod_lm$residuals^2)) #mse

[1] 8.414306

Se puede apreciar que el mse y msle es muy grande, el modelo no funciona bien.

plot(p1_train$count)
abline(mod_lm)

Se puede ver cómo un modelo lineal es demasiado sencillo para entender la complejidad de los datos y senzillamente predice la media del resultado.

Podremos coger el error 20000 (que representa el mse entre la media y cada valor) para poder decidir en los modelos posteriores si estan aprendiendo.

Validación sobre el test

preds_test_lm <- predict(mod_lm, test)
predsLM <- data.frame(pred = preds_test_lm, gt = test$count)
predsLM

Se puede apreciar cómo era de esperar, que la performance del modelo sobre el dataset de validación es mala.

Calculamos el MSE sobre el test

mean((predsLM$pred-predsLM$gt)^2) #mse

[1] 19119.78

Calculamos RMSLE (error cuadratico medio logaritmico) sobre el test

rmsle(predsLM$gt,predsLM$pred)

[1] NaN

4.3 glm

mod_glm <- glm(count~., train, family = poisson())
mean(mod_glm$residuals^2)

[1] 0.7673455

Se ha intentado analizar la predicción pero no se encuentra sentido en la misma. El mse de 0.79 parece demasiado bajo para ser correcto y el predict sobre el train da valores count lejanos al ground truth.

4.4 Random forest

set.seed(1234) #definimos una random seed para obtener los mismos resultados consistentemente

Definimos el modelo. Durante la realización de la práctica se ha hecho un grid search y se definen los parametros que devuelven mejor resultado.

rf <- randomForest(count~., data = train, mtry = 10, importance = TRUE, ntree = 50 )
rf

Call:
 randomForest(formula = count ~ ., data = train, mtry = 10, importance = TRUE,      ntree = 50) 
               Type of random forest: regression
                     Number of trees: 50
No. of variables tried at each split: 9

          Mean of squared residuals: 2346.799
                    % Var explained: 93.01

Visualizamos que predictores son más importantes.

importance(rf)

              %IncMSE IncNodePurity
year        74.922249      18359755
hour       163.338072     123392475
season      17.796134       7446561
holiday      8.081658        720623
workingday  87.259677      15505385
weather     19.430596       3216781
temp        47.551007      25933109
humidity    27.728905       7634230
windspeed    6.629868       2741114

varImpPlot(rf)

Podemos apreciar cómo la variable ‘hour’ es con diferencia la que mejor ayuda a predecir el numero de bicicletas alquiladas.

Analisis del model performance (mse) sobre el conjunto de train

mseDF <- data.frame(pred = rf$predicted, gt = rf$y)
mseDF

Podemos ver las diferencias entre el count predecido y el ground truth.

mean((mseDF$pred-mseDF$gt)^2) #mse

[1] 2346.799

Calculamos RMSLE (error cuadratico medio logaritmico)

rmsle(mseDF$gt,mseDF$pred)

[1] 0.3805507

Obtenemos un mse y rmsle coherentes, el modelo funciona.

Validacion del modelo sobre el test

pred_test_rf <- predict(rf, test)
mseDFtest <- data.frame(pred = pred_test_rf, gt = test$count)
mseDFtest

Podemos ver las diferencias entre el count predecido y el ground truth.

mean((mseDFtest$pred-mseDFtest$gt)^2) #mse

[1] 2133.669

Calculamos RMSLE (error cuadratico medio logaritmico)

rmsle(mseDFtest$gt,mseDFtest$pred)

[1] 0.3766893

El resultado del mse y el rmsle sobre el dataset de test es muy pareciodo a los resultados de train, de esta manera podemos concluir que nuestro modelo no tiene overfitting y es capaz de generalizar y tener buen rendimiento predictivo sobre datos que no ha visto con anterioridad.

4.4.1 Random Forest 2

Volvemos a hacer random forest con variables más correlacionadas con la columna a predecir.

trainImp <- train %>% select('hour','year','workingday','season','humidity', 'count')

Entrenamos random forest

rf2 <- randomForest(count~., data = trainImp, mtry = 10, importance = TRUE, ntree = 50 )
rf2

Call:
 randomForest(formula = count ~ ., data = trainImp, mtry = 10,      importance = TRUE, ntree = 50) 
               Type of random forest: regression
                     Number of trees: 50
No. of variables tried at each split: 5

          Mean of squared residuals: 4280
                    % Var explained: 87.26

Comparación predicted vs ground truth.

mseDF2 <- data.frame(pred = rf2$predicted, gt = rf2$y)
mseDF2

Calculamos error

mean((mseDF2$pred-mseDF2$gt)^2) #mse

[1] 4280

Fucniona peor cogiendo solo las variables más importantes (mse 4280 > mse 2346). Nos quedaremos con el modelo ‘RF’. Al tener peor performance sobre el train, no se haran las pruebas sobre el dataset de test ya que se asume que el modelo no rendira mejor que en el caso anterior.

5. Test

Realizaremos una predicción sobre el dataset de test con el mejor modelo, en nuestro caso el RF para la entrega final.

5.1 Realizar las mismas transformaciones en el test que en el train

p1_test <- p1_test %>% mutate(workingday = as.factor(workingday),
                                   holiday = as.factor(holiday),
                                   weather = as.factor(weather),
                                   season = as.factor(season)) %>% mutate(across(where(is.numeric), scale))

5.2 Realizar la prediccion

Se realiza la prediccion sobre los datos de test con el modelo RandomForest.

predsFin <- predict(rf, p1_test)
predsFinDF <- data.frame(id = p1_test_fin$id, prediction_test = predsFin)
predsFinDF

Se genera el DF para visualizar predicciones.

5.3 Exportar resultados

Se exportan los resultados con el id y el resultado siguiendo las instrucciones del enunciado.

write.table(x = predsFinDF, file = 'p1.txt', sep = ',', row.names = FALSE, col.names = TRUE)

3. Conclusiones

El algoritmo de regresion random forest es el que mejor predice el numero de bicicletas que se alquilaran un dia, con un mse de 2000 sobre el conjunto de train. Podemos ver analizando los resultados que el modelo no tiene mucho overfitting y generaliza bien, esto es asi en parte al haber creado un arbol no demasiado grande, evitando así la caída de valores en hojas.

Teniendo en cuenta que el mse entre la media y todos los valores del dataset era de 20.000 y con el mejor modelo hemos obtenido un mse de 2000, podemos concluir que nuestro modelo ha aprendido satisfactoriamente.

El rmsle definitivo es 0.3695941.

Cómo algoritmos de regresion se podrian intentar también el xgboost o el catboost pero perderiamos entendimiento de lo que esta haciendo el modelo, asi que se ha obviado para este trabajo.